Fine Structure of Oscillons in the 
Spherically Symmetric Klein-Gordon Model 

Ethan P. Honda 

Applied Research Laboratories, The University of Texas at Austin, Austin, TX 78758 

Matthew W. Choptuik 
CIAR Cosmology & Gravity Program 
Department of Physics and Astronomy, 
University of British Columbia, 
Vancouver BC, V6T IZl 
Center for Relativity, 
University of Texas at Austin, TX 78712-1081 USA 

We present results from a study of the fine structure of oscillon dynamics in the 3+1 spherically 
symmetric Klein-Gordon model with a symmetric double-well potential. We show that in addition 
to the previously understood longevity of oscillons, there exists a resonant (and critical) behavior 
which exhibits a time-scaling law. The mode structure of the critical solutions is examined, and we 
also show that the upper-bound to oscillon formation (in ro space) is either non-existent or higher 
than previously believed. Our results are generated using a novel technique for implementing non- 
reflecting boundary conditions in the finite difference solution of wave equations. The method uses a 
coordinate transformation which blue-shifts and "freezes" outgoing radiation. The frozen radiation 
is then annihilated via dissipation explicitly added to the finite-difference scheme, with very little 
reflection into the interior of the computational domain. 

PACS numbers: 04.25.Dm,04.40.Nr,11.10.Lm,11.27.-f d,64.40.Ht,98.80.Cq 



I. INTRODUCTION 



There is a long history in physics and mathematics of trying to find new non-trivial solutions to nonlinear wave 
equations. The literature on the subject goes back at least as far as 1845 when J. Scott Russell published a paper 
about a surface wave he witnessed travel for almost two miles in a shallow water channel (the first scientifically 
reported soliton)[l]. Since then there has been much effort directed towards understanding stable localized solutions 
to nonlinear wave equations: the classical kink-soliton, topological defects (monopoles, cosmic strings, and domain 
walls) [2], and nontopological defects (such as Q-balls)[3, 4], are but a few examples. However, localized but unstable 
solutions have been discussed much less frequently and in this paper we focus our attention on one such solution, the 
oscillon. 

The definition of oscillon varies slightly depending on context, but here we refer to localized, time-dependent, 
unstable, spherically symmetric solutions to the nonlinear Klein-Gordon equation. Although oscillons are unstable, 
their lifetimes are long compared to a dynamical time. Oscillons were originally discovered by Bogolubsky and 
Makhan'kov [5, 6] (who called them "pulsons"), and were later studied in more detail by Gleiser [7] and by Copeland 
et al [8]. 

Oscillons can be formed via the collapse of field configurations (initial data) that interpolate between two vacuum 
states (0+ and 0_) of a symmetric double well potential (SDWP)^. In spherical symmetry, such a configuration is a 
bubble, where the interpolating region is the bubble "wall" that separates the two vacuum states at some characteristic 
radius (where in this work we always use as the large-r vacuum state) . An oscillon formed this way typically has 
three distinct stages in its evolution. First, immediately following the bubble collapse a large percentage of its energy 
is radiated away. As will be discussed below, this can happen either through localized field oscillations, or through 
bounces reminiscent of the 1-1-1 dimensional kink-antikink, (KK), scattering [9]. After the initial radiative phase, 
the solution settles into the oscillon stage. Here the field is localized with a shape roughly that of an origin-centered 
gaussian, with the field value asymptotically approaching the large-r vacuum state. Due to the asymmetry of the 
potential about either minimum, the field oscillates about 4>- (the large-r vacuum) such that the time-averaged value 



^ Asymmetric double well potentials can also produce oscillons, but here we consider only the SDWP. 
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of the field lies between the two vacua, ic. 4)- < {4>) < (t>+ where (•••)= ^'"Vo ' ' ' ^* 

as in [8]. For typical initial 

field configurations the energy of the oscillon is slowly radiated away, approaching a specific "plateau" value. In the 
third and final stage of evolution, the oscillon stops oscillating and disperses, radiating away its remaining energy. 

Much of the original excitement about oscillons arose from the fact that their long lifetimes could potentially alter 
the dynamics of a cosmological phase transition. However, since oscillons are unstable, their ability to affect such 
a phase transition depends crucially on their lifetimes. Previous studies by Copeland et al [8], used dynamical grid 
methods to study oscillon dynamics and treated the initial radius, ro, of the bubble (shell of radiation), as a free 
parameter. These studies showed not only that oscillon lifetimes can be comparable to the age of the universe (at 
the GUT scale), but that oscillons are formed from a wide range of initial bubble radii. However, the computational 
demands of the dynamical grid methods employed in [8] prevented a detailed study of the ro parameter space. 

A key problem in the accurate, long-time simulation of oscillons is the treatment of boundary conditions at the outer 
edge, r = rmax. of the computational domain. It is standard practice in the computational solution of nonlinear field 
equations to use finite difference techniques applied to functions defined on a lattice of gridpoints. If a static, finite- 
sized domain is used; i.e. if r^ax is fixed, then one needs to employ a method which minimizes the amount of radiation 
(energy) which is artificially reflected at r = rmax- With massless scalar fields, and in spherical symmetry, this can 
be done quite easily simply by imposing a discrete version of an "outgoing radiation" , or Sommerfeld, condition. 
However, for the case of massive scalar fields, or more generally, for fields with non-trivial dispersion relations, the 
Sommerfeld condition is only approximate, and its use generically results in significant reflection at r = rmax, and 
subsequent contamination of the interior solution. 

A surefire fix for the outer-boundary problem is to use a dynamically growing grid (as in [8]), so that rmax = rmax(i), 
and lattice points are continuously added to extend the computational domain as needed. Alternatively, compactified 
coordinates, or coordinates which propagate outwards faster than any characteristic speed in the problem can be used, 
but in these cases, new gridpoints also need to be continuously added to the mesh in order to maintain adequate 
resolution of solution features. These methods are somewhat more efficient than the use of a static mesh with rmax 
chosen so that no signals reach the outer boundary during the integration period of interest, T. However, for long-lived 
solutions, the mesh soon becomes quite large, and the computation time tends to be proportional to T^. 

Recently, Gleiser and Sornborger [10] introduced an adiabatic damping method which adds an explicit damping term 
to the equations of motion, and which has been shown to absorb outgoing massive radiation extremely well in ID 
(spherical) and 2D (cylindrical) simulations. Here we present an alternate approach for dealing with outgoing massive 
scalar fields which is quite general and quite different from previously used methods of which we are aware. The 
technique involves the use of a specially chosen coordinate system that "freezes" and blue-shifts outgoing radiation 
in a relatively thin layer located away from the central region where the dynamics of principal interest unfold. The 
addition of a standard type of finite-difference dissipation [11] then "quenches" the blue-shifted, frozen radiation, 
and very little energy is reflected back into the interior region. This approach, like that described in [10], has the 
advantage that a static and uniform finite-difference mesh can be used, so that computational time scales linearly 
with the integration period, T. 

Our new techniqTie was crTicially important to our discovery and detailed study of fine striicture in a well-known 
(and still much studied) nonlinear system. Specifically, we have found strong evidence for a family of resonant oscillon 
solutions in the SDWP model. Each of these solutions appears to possess a single unstable mode in perturbation 
theory, and by tuning the family parameter, ro, in the vicinity of a specific resonance, we can tune away that mode, 
producing oscillons which live longer and longer as we tune closer and closer to the precise resonant value, Tq. This 
leads to a view of oscillons as being analogous to the Type I critical solutions which have been discovered in the 
context of gravitational collapse [12], and as in that case, we find compelling evidence for power-law scaling of the 
oscillon lifetime, r: 

r ~ Cr\ra - r*\'''- (1) 

where Cr is an overall scale factor set by the particular resonance, and 7^ is a resonance- dependent exponent which is 
presumably the reciprocal Lyapounov exponent associated with the resonance's single unstable mode. 

In addition, contrary to previous claims [8, 10], we see no hard evidence for an upper bound on ro, beyond which 
oscillons are no longer generated via collapse of gaussian data. In particular we find strong evidence for resonances 
for ro ^ 6.5, well above the limit ro ~ 4.2 quoted in [8, 10]. Moreover, we relate the existence of these "large-ro" 
resonances to the "bouncing" behaviour observed in the 1-1-1 kink-antikink study of Campbell et al [9]. 

The remainder of the paper is organized as follows: In Section 2 we introduce a new coordinate system in which 
to solve non-linear wave equations using finite differences. We examine the conformal structure induced by our new 
coordinates, as well as the characteristics of the resulting wave equation. In Section 3 we discuss the new properties 
of oscillons which were discovered during our study. In particular, we observe resonances in the parameter space 
which obey a time-scaling law, and we construct a sample resonant solution via a non-radiative ansatz (Sections 3a 
and 3b, respectively). Finally, in Section 3c we discuss oscillons and resonant solutions found outside the bounds of 
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the parameter space previously explored. Section 4 summarizes our results and is followed by two appendices which 
discuss the details of the finite difference equations (appendix A) and the testing of the code (appendix B). 

II. THE KLEIN-GORDON EQUATION IN MIB COORDINATES 

We are interested in the self-interacting scalar field theory described by the (3+l)-dimensional action 

5[</.] = j d'x^\ (-\g'"'V^ct^V,ct^ - Vict^)^ (2) 
where we take V{(j}) to be a symmetric double well potential (SDWP) ^, 

Vs{d,) = \{<f-lf (3) 

and g^v to be the metric of flat spacetime in spherical symmetry, written in standard spherical polar coordinates 

{t,f,e, If): 



ds 



p:2 



= -dp + df^ + ^ gjn^ e #2 j (4) 



We now introduce a new radial coordinate, r, which interpolates between the old radial coordinate, r, at small r and 
an outgoing null coordinate at large f. Specifically, we take 

t = t (5) 

f = r + f{r)t (6) 

= 9 (7) 

^ = ip (8) 

where /(r) is a monotonically increasing function which smoothly interpolates between « and « 1 at some character- 
istic cutoff, Tc, so that /(r) — > for r <C Tc, and /(r) ^ 1 for r rc- We call (t,r) monotonically increasingly boosted 
(MIB) coordinates. The MIB system reduces to the original spherical coordinates, (t,r), for r -^r^, but as discussed 
below, in the r > rc region, both outgoing and ingoing (from r ^ rc) radiation tends to be "frozen" in the transition 
layer, r « rc- Furthermore, since the outgoing radiation is blue-shifted as it propagates into the transition region, 
r rc, application of standard finite-difference dissipation operators can then quench it with minimal reflection back 
into the interior of the computational domain. 

In general, MIB coordinates will not cover all of the {i,f) half-plane. However, given that /(r) is monotonically 
increasing, the determinant of the Jacobian of the transformation is non-zero for all t such that t > — max|/'(r)|. 
Thus, for this range of t, the transformation to and from the standard spherical coordinate system is well-defined, 
and though a coordinate singularity inevitably forms as t — oo (past timelike infinity), this has no effect on the 
forward evolution of initial data given at t = 0. 

We also note that in order that our MIB coordinates be regular at r = (so that there is no conical singularity at 
the origin), we must also demand that /(O) = 0. 

Our coordinate choice results in the following spherically symmetric, 3-1-1, or ADM [13] form: 

ds^ = {-a^ + a^(3^) dt" + 2a^f3dtdr + a^dr^ + r^fe^ (^^2 ^ gjj^2 g ^^2^ 

where 



a(t,r) = l + /'(r)f &(t, r) = l + /(r)- 

r 

f{r) 



a{t,r) = 1 I3{t,r) 



(10) 



1 + nr)t 



^ This is identical to using V(^) ~ '^{^^ ^ ] introducing dimensionless Vciriables r = fm, t = im, and x = '^'t'- 
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(In the nomenclature of the ADM formahsm. a is the lajtse junction^ while /? is the radial component of the shift 
vector.) In the work which follows, we have adopted the following specific form for /(r): 

/(r) = (1 + tanh((r - r,)/6)) /2 + e, (11) 

where 

e = - (1 + tanh(rc/5)) /2 (12) 

is chosen to satisfy the regularity condition at r = 0. 

It is now instructive to consider the conformal structure of the MIB hypersurfaces. This is done by applying 
equations (8) to the standard conformal compactification on Minkowski space, i ± f = tan((r±i?) /2) (where T 
and R are the axes in the conformal diagram, see [14] or [15]), and then plotting curves of constant r and t. The 
constant-t hypersurfaces are everywhere spacelike and all reach spatial infinity, i°. Although constant-r surfaces for 
r > Tc appear at first glance to be null, a closer look (see insets of Fig. 1) reveals that they are indeed everywhere 
timelike and do not ever reach future null infinity, JZ"*". 

The equation of motion for the scalar field which results from the action (2) is 

-^d^ (VW\9'"'d.A = </. - 1) (13) 
V\9\ ^ ' 

which with (9), (10), 11 = a (9(0 — /39r0) /a, and $ = 9r0 give 



n 



a 

I 



^ (r%- {U + /3n))' - 2^n - aa0 (0^ - l) (14) 

(^n + (15) 
"n + /3$ (16) 



where ' = dt and ' = d^. These equations are familiar from the ADM formalism as applied to the spherically 
symmetric Klein-Gordon field coupled to the general relativistic gravitational field [16] . However, in the current case, 
instead of a dynamically evolving metric functions, the metric components a, h, a, and (3 are a priori fixed functions 
of [t, r) that resulted from a coordinate transformation of flat spacetime. 

Clearly, characteristic speeds for the masslcss Klein-Gordon field (V(0) = 0) bound the inward or outward speed 
(group velocities) of any radiation in a self-interacting field (V^(0) 7^ 0). Characteristic analysis of the masslcss 
Klein-Gordon equation with metric (9) yields local propagation speeds 

A± = -/3±^, (17) 

where A+ and A_ are the outgoing and ingoing characteristic speeds, respectively [16, 17]. For r -C rc , the propagation 
of scalar radiation in (i, r) or {i,f) coordinates is essentially identical. However, as illustrated in Fig. 3, and as can 
be deduced from equations (17) and (10), for r « rc both the ingoing and the outgoing characteristic velocities go to 
zero as t — > cxD (as the inverse power of t). Thus, any radiation incident on this region will efi'ectively be trapped, 
or "frozen in". It is this property of the MIB system which enables the effective implementation of non-reflecting 
boundary conditions. As discussed further in Appendix A, an additional key ingredient is the application of Kreiss- 
Oliger-style dissipation [11] to the difference equations. This dissipation efficiently quenches the trapped outgoing 
radiation, which as mentioned above tends to be blue-shifted to the lattice scale on a dynamical time-scale. 

Finally, we note that, as is evident from Fig. 3, the "absorbing layer" in the MIB system (i.e. the region in which 
the characteristic speeds are « 0), expands both outward and inward as t increases. This means that for fixed rc, 
the absorbing layer will eventually encroach on the interior region r <§; and ruin the calculation. However, the 
rate at which the layer expands is roughly logarithmic in t, so, in practice, this fact should not significantly impact 
the viability of the method. For arbitrarily large final integration times, T, computational cost will scale as TlnT. 
However, the calculations described here all used the same values of Tc and Tmax, so that for all practical purposes, 
the computational cost is linear in the integration time. 



III. THE RESONANT STRUCTURE OF OSCILLONS 



Copeland, et al, showed quite clearly that oscillons formed for a wide range of initial bubble radii, tq. They even 
caught a glimpse of the fine structure in the model — which in large part motivated this study — but they did not 
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explore this fine structure of the parameter space in detail. With the eSiciency of our new code, we have been able 
to explore parameter space much more thoroughly, which in turn has yielded additional insights into the dynamical 
nature of oscillons. 

Following [7] we use a gaussian profile for initial data where the field at the core and outer boundary values are set 
to the vacuum values, (j){t, 0) = </)c = 1 and (j){t, oo) = = — 1 respectively, and the field interpolates between them 
at a characteristic radius, ro- 

(f>{0, r) = 0o + (<^e - M exp {-r^/rl) . (18) 

Keeping (j)c and <po constant, but varying ro, we have a one parameter family of solutions to explore. Figure 4 shows 
the behavior of oscillon lifetime as a hmction of ro in the range 2.0 < ro < 5.0. We discuss three main findings that 
are distinct from previous work: the existence of resonances and their time scaling properties, the mode structure of 
the resonant solutions, and the existence of oscillons outside the parameter-space region 2 < ro < 5. 

A. Resonances & Time Scaling 

In contrast to Fig. 7 of Copeland et al [8], the most obvious new feature seen in Fig. 4 is the appearance of the 125 
resonances which rise above the overall lifetime profile. These resonances (also seen in Fig. 5) become visible only 
after carefully resolving the parameter space. Upon fine-tuning ro to about 1 part in ~ 10^^ we noticed interesting 
bifurcate behavior about the resonances (Figure 6, top). The field oscillates with a period T « 4.6 (for all oscillons) 
so the individual oscillations cannot be seen in the plot, but it is the lower-frequency modulation that is of interest 
here^. The top figure shows the envelope of (/)(t, 0) on both sides of a resonance (dotted and solid curves). We see that 
the largo period modulation that exists for all typical oscillons disappears late in the lifetime of the oscillon as ro is 
brought closer to a resonant value, r^. On one side of the modulation returns before the oscillon disperses (refered 
to as supercritical and shown with the solid curve), while on the other side of rj the modulation does not return 
and the the oscillon simply disperses (refered to as subcritical and shown with the dotted curve). For resonances 
where ro ^ 2.84, the subcritcal solutions appear on the ro < rp side of the resonance and the supercritical solutions 
appear on the ro > r^ side of the resonance. The opposite is true for resonances where r^ > 2.84, ie. the subcritcal 
solutions appear on the ro > r^* side of the resonance and the supercritical solutions appear on the ro < rg side of 
the resonance. This bifurcate behavior does not manifest itself until ro is quite close to r^. In practice then, to locate 
a resonant value, rj, we first maximize the oscillon lifetime using a three point extremization routine {golden section 
search with bracketing interval of ~0.62, [18]) until we have computed an interval whose end-points exhibit the two 
distinct behaviours just described. Once a resonance has been thus bracketed, we switch to a standard bisection 
search, and subsequently locate the resonance to close to machine precision. Although we can see from Fig. 6 that 
the modulation is directly linked to the resonant solution, it is not obvious why this is so. However, if we look at the 
relationship between the modulation in the field (top) to the power radiated by the oscillon (bottom), we see that 
they are clearly synchronized. 

The behaviour of these resonant solutions may not be surprising to those familiar with the 1+1 KK scattering 
studied using the same model [9] . Campbell et al showed that after the "prompt radiation" phase — the initial release 
of radiation upon collision of a kink and anti-kink — the remaining radiation was emitted from the decay of what they 
referred to as "shape" oscillations. The "shape modes" were driven by the contribution to the field "on top" of the 
K and K soliton solutions. Since the exact closed-form solution for the ideal non-radiative KK interaction is not 
known, initial data aimed at generating such an interactionm is generally only approximate, and the "surplus" (or 
deficit) field is responsible for exciting the shape modes. The energy stored in the shape modes slowly decays away 
as the kink and antikink interact and eventually the solution disperses. 

In our case, wc believe the large period modulation represents the excitation of a similar "shape mode" superimposed 
on a periodic, non-radiative, localized oscillating solution. On either side of a resonance in the ro parameter space, the 
solution is on the threshold of having one more shape mode oscillation. If this is the case, then, as we tune ro ^ Tq, 
we are, in effect, tuning away the single unstable shape mode, and thus should expect that the oscillon lifetime obey 
a scaling law such as that seen in Type I solutions in critical gravitational collapse [12]. Figure 7 shows a plot of 
oscillon lifetime versus In |ro — rgj (for the ro = 2.335 • • • resonance), and we can see quite clearly that there is a 
scaling law, T ~ 7ln |ro — rgl, for the lifetime of the solution as measured on either side of the resonance. We denote 
7+ for the scaling exponent on the ro > ro side, and 7_ for the scaling exponent on the ro < ro side. Although we 



In dimcnsionful coordinates, f and t, the period would be T = 4.6m ^ . In general, to recover proper dimensions, lengths and times are 
multiplied by and energies by \m~^. 
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observe lifetime scaling for each resonance, the scaling exponent per se varies from resonance to resonance; a plot of 
the scaling exponents, 7+ and 7_, versus the critical initial bubble radius can be seen in Fig. 8. For all resonances 
we find 7_|_ » 7_ . 

Finally we note that, by analogy with the case of Type I critical gravitational collapse, we expect that the scaling 
exponents, 7, are simply the reciprocal Lyapounov exponents associated with each resonance's single unstable mode. 
In addition we note that, for any resonance, if we were able to infinitely fine-tune ro to Tq, we would expect the 
oscillon lifetime to go to infinity. 



B. Mode Structure 



Assuming that periodic, non-radiative solutions to equation (13) exist, we should be able to construct them by 
inserting an ansatz of the form 

00 

(l){t, r) = Mr) + ^"('') ("^*) (^^) 

n=l 

in the equations of motion and solving the resulting system of ordinary differential equations obtained from matching 
cos {nuit) terms: 



+i(<^o-l)E('Am)' 



4 

m,p,q 



(20) 



(rX)'/r-' = (3 (00 + 

+ l{(l)0-l)Y.(pp(pq{5n,±p±q) I2l) 
p,q ^ I 

+ 2 Yli <t'm<t)n<t>q {5n,±m±p±q) ■ 
m,p,q 

Equations (20) and (21) can also be obtained by inserting ansatz (19) into the action and varying with respect to 
the (t)n [19]. This set of ODEs can be solved by "shooting", where the quantities 0n(O) are the shooting parameters. 
Unfortunately, we were unable to construct a method that self-consistently computed w; the best we could achieve 
was to solve equations (20) and (21) for a given u which we measured from the PDE solution. 

For ease of comparison of the results obtained from the periodic ansatz with those generated via solution of the 
PDEs, we Fourier decomposed the PDE results. This was done by taking the solution during the interval of time when 
the large period modulation disappears (1200 < t < 1800 for the oscillon in Fig. 6, for example) and constructing 
FFTs of (j) at each gridpoint, r^. Specifically, at each r^, the amplitude of each Fourier mode was obtained from a 
FFT which used a time series, (t)(t", ?';), n = 1,2, • • • 4096 with f^^^ —t" = At = constant. Keeping only the first five 
modes in the expansion (19), we compare the Fourier decomposed PDE data with the shooting solution (see Fig. 9). 
It should be noted that although the value for ui was determined from the PDE solution, the shooting algorithm still 
involved a five-dimensional search for the the shooting parameters, (/)„(0), n = 0, • • • , 4. The close correspondence of 
the curves shown in Fig 9 strongly suggests that the resonant solutions (ie. in the limit as ro ^ Tq) observed in the 
PDE calculations are indeed consistent with the periodic, non-radiative oscillon ansatz (19). 

By examining the three most dominant components of the power spectrum of d>{t, 0), Fig. 10, we can see that during 
the "no-modulation" epoch, the amplitude of each Fourier mode becomes constant. Although the specific plot is for 
the core amplitude, r = 0, we note that this behavior holds for all r. Again, this is consistent with the view that as we 
tune ro to rj, the oscillon phase of the evolution is better and better described by a one- mode unstable, "intermediate 
attractor" . As discussed previously, this is precisely reminiscent of the Type I critical phenomena studied in critical 
gravitational collapse, particularly the collapse of a real, massive scalar field as studied by Brady et al [12], where the 
intermediate attractors are unstable, periodic, "oscillon stars" discovered earlier by Seidel and Suen [20]. 



C. (Bounce) Windows to more Oscillons 

Lastly, we consider the existence of oscillons generated by gaussian initial data with ro ^ 5. The oscillons explored 
by Copeland et al, were restricted to the parameter-space region, 2 < ro ^ 5, and in fact it was concluded that there 
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was an upper bound, tq ^ 4.2, beyond which evolution of gaussian data would not result in an oscillon phase [10]. 
However, we have found that oscillons can form for tq 5, and that they do so by a rather interesting mechanism. 

Again, from the 1+1 dimensional KK scattering studies of Campbell et al, it is well known that a kink and antikink 
in interaction often "bounce" many times before either dispersing or falling into an (unstable) bound state. Here, 
a bounce occurs when the kink and antikink reflect off one another, stop after propagating a short distance, then 
recollapse. 

We find that such behaviour occurs in the (3+1) dimensional case as well, but now the unstable bound state is an 
oscillon. For larger ro, instead of remaining within r < 2.5 after reflection through r = (as occurs for 2 < ro ^5), 
the bubble wall travels out to larger r (typically 3 r ^ 6), stops, then recoUapses, shedding away large amounts of 
energy in the process (see Fig. 11). Thus in this system, as with the 1+1 KK model, there are regions of parameter 
space which constitute "bounce windows" . Within such regions, the bounces allow the bubble to radiate away large 
amounts of energy. The bubble then recoUapses, effectively producing a new initial configuration (albeit with a 
different shape) with a smaller effective tq. Within these "windows" both oscillons and resonances (similar to those 
observed for 2 < ro 4.6) can be observed (inset of Fig. 12). 

IV. CONCLUSIONS 

Using a new technique for implementing non-reflecting boundary conditions for finite-differenced evolutions of non- 
linear wave equations, we have conducted an extensive parameter space survey of bubble dynamics described by a 
spherically-symmetric Klein-Gordon field with a symmetric double- well potential. We have found that the parameter 
space of the model exhibits resonances, wherein the lifetimes of the intermediate-phase "oscillons" diverge as one 
approaches a resonance. We have conjectured that these resonances are single-mode unstable solutions, analogous to 
Type I solutions in critical gravitational collapse, and have presented evidence that their lifetimes satisfy the type of 
scaling law which is to be expected if this is so. 

In addition, we have independently computed resonant solutions starting from an ansatz of periodicity, and have 
demonstrated good agreement between the solutions thereby computed, and those generated via finite-difference 
solution of the PDEs. Finally, we have showed that oscillons can form from bubbles with energies higher than 
had previously been assumed, through a mechanism analogous to the bounce windows found in the 1 + 1 case of 
kink-antikink scattering. 

We note that the use of MIB or related coordinates, in conjunction with finite-difference dissipation techniques, 
should result in a generally-applicable strategy for formulating non-refiecting boundary conditions for finite-difference 
solution of wave equations. The method has already been used in the study of axisymmetric oscillon collisions [21], 
and attempts are underway to use similar techniques in the context of 3-D numerical relativity and 2-D and 3-D ocean 
acoustics. 
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APPENDIX A: FINITE DIFFERENCE EQUATIONS 

Equations (14,15,16) are solved using two-level second order (in both space and time) finite difference approximations 
on a static uniform spatial mesh: 

ri = {i-l)Ar, i = l,2,---7 (Al) 
where / is the total number of mesh points 

= . (A2) 

The scale of discretization is set by Ar and At = AAr, where we fixed the Courant factor. A, to 0.5 as we changed 
the base discretization. 

Using the operators from Table I, drf = a, dr = nr"'~^drn, and rb = f, the difference equations applied in the 
interior of the mesh, i = 2, 3, • • • 7 — 1, are 

AfU2 = S/iit [aAf.3 (f2 (-$ + /3n))]" 

fb (A3) 
-2/x J -n - aa(^ ((/)2 - 1) J , 

A,^$r = MtA, (-n + /3$)", (A4) 

A^0r = Mt(-n + /3$)". (A5) 

These equations arc solved using an iterative scheme and explicit dissipation of the type advocated by Krciss and 
Oliger [11]. The dissipative term, incorporated in the operator A^ , is essentially a fourth spatial derivative multiplied 
by (Ar)^ so that the truncation error of the difference scheme remains 0{Ar^,At'^). The temporal difference operator, 
Af , is used as an approximation to di everywhere in the interior of the computational domain, except for next- to- 
extremal points, where Aj is used because the grid values /^'^j '^^ fi-2 ^^i* defined. 
At the inner boundary, r = 0, we use O(Ar^) forward spatial differences to evolve 11 

Afn-^nj =0. (A6) 

whereas is fixed by regularity: 

<E>? = 0. (A7) 
To update (j), we use a discrete versions of the equation for (j) which follows from the definition of 11: 

At<P^ ^J^t(—+|3^\ i=l,2,---I (AS) 



a 

At the outer boundary, r = rmax, our specific choice of boundary conditions and discretizations thereof have little 
impact; due to the use of MIB coordinates and Kreiss-Oliger dissipation, almost none of the outgoing scalar field reaches 
the outer edge of the computational domain. Nevertheless, we imposed discrete versions of the usual Sommerfeld 
conditions for a massless scalar field on 11 and $: 

Atn?+MtfA^n+5^ =0, (A9) 



At^^ + ^W A^$ + ^ ) =0, (AlO) 
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APPENDIX B: TESTING THE MIB CODE 



One might think that "freezing" outgoing radiation on a static uniform mesh would lead to a "bunching-up" of the 

wave-train from the oscillating source, which would then result in a loss of resolution, numerical instabilities, and an 
eventual breakdown of the code. However, this turns out not to be the case; all outgoing radiation is "frozen" around 
r ~ Tc, but the steep gradients which subsequently form in this region are efficiently and stably annihilated by the 
dissipation which is explicitly added to the difference scheme. 

In fact, there is a loss of resolution and second order convergence for r ~ rc, but this does not affect the stability 
or convergence of the solution for r <C Vc- Figure (14) shows a convergence test for the field for r < rc/2 over 
roughly six crossing times. Since we are solving equation (13) in flat spacetime, it is very simple to monitor energy 
conservation. The spacetime admits a timelike Killing vector, t" , so we have a conserved current, = t'^T^jy. We 
monitor the flux of through a surface constructed from two adjacent spacclike hypersurfaces for r < Tc (with 
normals = (±1, 0, 0, 0)), and an "endcap" aXr = ro (with normal = (0, a~^, 0, 0)). To obtain the the conserved 
energy at a time, tf, the energy contained in the bubble, 

ro 2 2 

i^bubbie = 4X252 +^ + v{ct>)\ dr, (Bl) 







2a2 

(where the integrand is evaluated at time tf) is added to the total radiated energy, 



b 



^JrH^^dt (B2) 



(where the integrand is evaluated at r = ro)- The sum, E'totai = -^bubble + E^ad, remains conserved to within a few 
tenths of a percenf* through a quarter million iterations (see Fig. (15)). 

Although monitoring energy conservation is a very important test, it says little about whether there is reflection 
of the field off of the outer boundary, r = Tmax, or the region r w Tc. To check the efficacy of our technique for 
implementing non-reflecting boundary conditions, we compare the MIB results to those obtained with two other 
numerical schemes. The flrst alternate method involves evolution of equation (13) in {t,f) coordinates on a grid with 
fmax sufficiently large that radiation never reaches the outer boundary (large-grid solutions). For a given discretization 
scale, results from this approach serve as near-ideal reference solutions, since the solution is guaranteed to be free 
of contamination from reflection off the outer boundary. The second method involves evolution on a grid with the 
same rmax adopted in the MIB calculation, but with discrete versions of massless Sommerfeld (outgoing radiation) 
conditions applied at r = rmax- We refer to the results thus generated as OBC solutions, and since we know that 
these solutions do have error resulting from reflections from r = r,nax, they demonstrate what can go wrong when a 
solution is contaminated by reflected radiation. Treating the large-grid solution as ideal. Fig. (16) compares typical 
logio 11?^ ~ (/"ideailb for the MIB and OBC solutions. There is a steep increase in the OBC solution error (three orders 
of magnitude) around t = 125, which is at roughly two crossing times. This implies that some radiation emitted 
from the initial collapse reached the outer boundary and reflected back into the region r < ro- There is no such 
behavior found in any MIB solutions. Lastly, for a more direct look at the held itself, we can see <j>{t,0) for large-grid 
(triangles), MIB (solid curves), and OBC (dashed curves) solutions in Fig. (17). Initially, both the MIB and OBC 
solutions agree with the large-grid solution extremely well. However, after two crossing times the OBC solution starts 
to substantially diverge from the ideal solution, while the MIB results remain in very good agreement with the ideal 
calculations. 

In summary, the MIB solution conserves energy, converges quadratically in the mesh spacing (as expected), and 

produces results which are equivalent — at the level of truncation error to large-grid reference solutions. At the same 
time, the MIB approach is considerably more computationally efficient than dynamical- or large-grid techniques. 
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Operator Definition Expansion 

Ajif (-3/r+4/i;i~/r+2)/2Ar a./|;+0(Ar2) 

A^/r (3/r - 4/r_i + /r_2) /2Ar a./|" + 0(Ar2) 

A./r (/r+1 - /r-i) /2Ar + O (Ar^) 

At/r (/r+' - /r) /At + o (At^) 

A*^/r (/r+i - /r) /At+ + o {/^e) 

edis[&fi + fT-2 + fi+2- 

4(/r_i+/r+i)]/16At 

^Hfl (/r+^+/r)/2 /|;+^+o(At^) 



TABLE I: Two-Level Finite Difference Operators. Here we have adopted a standard finite-difTerence notation: /" = f{{n — 
l)At, (i - l)Ar) 
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FIG. 1: Conformal diagram showing surfaces of constant r (dotted lines) and lines of constant t (dashed lines). Lines of constant 
t look exactly like the constant-t hypersurfaces of Minkowski space, whereas the lines of constant r behave much differently. 
For r > Tc, it appears as if the constant-r surfaces are null. However, as insets A, B, and C show, the constant-r lines do not 
become null (do not intersect future null infinity), but rather are everywhere timelike. 
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FIG. 2: Plot of the characteristic velocities as a function of the MIB coordinates, r' and t' in units where Vc is set to unity. A+ 
and A_ are the outgoing and ingoing characteristic speeds respectively. S is taken to be (5 — » S/rc « 0.0893 (corresponding to 

the system used in this article). Characteristic velocities arc plotted for times t' = 0, 1, 10, and 100 {t' = 100 is larger than 
the lifetime of the longest lived solution studied in this work). As r ^ 1, we have X±{t',r') = 0{l/t'). 
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FIG. 3: Plot of characteristic speeds, A±(r', 100), where r' and t' are radial MIB coordinates in units where Tc is set to unity. 
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FIG. 4: Plot of oscillon lifetime versus initial bubble radius for 2.0 < ro < 5.0. Each of the 125 resonances is resolved to one 
part in lO". 
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FIG. 5: Plot of oscillon lifetime versus initial bubble radius for 2.27 < rinitiai < 2.29. The three resonances shown occur 
at ~ 2.2805, rl ~ 2.2838, and r| ~ 2.2876. Each resonance separates the parameter space into regions with n and n + 1 
modulations; the x's correspond to oscillons with no modulations, the triangles to oscillons with one modulation, the squares 
to two modulations, and the circles to three modulations. 
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FIG. 6: Top plot shows the envelope of 4i{t,0) for rgiAro displaying bifurcate behavior around the rg « 2.335 resonance 
(Aro ~ 10""); the solid curve is the envelope barely above resonance (15 modulations) while the dotted line is the envelope 
barely below resonance (14 modulations). Bottom plot shows the energy radiated as a function of time through the surface 
containing the oscillon as defined in Appendix B. The increases in the energy radiated are synchronized with the modulation 
in the field. 
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FIG. 7: Plot of time scaling, T versus In |ro — rgl for the ro ~ 2.335 resonance. The top line (triangles) displays the scaling 
behavior for supercritical evolutions, ro > Tq, while the bottom line (x's) shows the scaling for subcritical calculations, ro < ro. 
The exponents (measured by the slopes of the lines) are both approximately equal to 7 = 30. The offset in the two curves 
represents the time spent by supercritical oscillons in executing the final modulation shown in Fig. 6. 
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FIG. 8: Plot of critical exponents for each resonance. There are two values of 7 for each resonance. The top plot displays the 
7+ vs. tq while the lower plot displays 7_ vs. rj. The uncertainties are estimated from running the entire parameter space 
surveys at two resolutions, A*' = A*',- = 1449 and TV' = N'^ — 1025 and estimating the error, A7 — \ jn — 7jv'|- 
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FIG. 9: Critical solution 4>n{r) (for n = 0, 1, 2, 3, 4) obtained from the Fourier-decomposed PDE data (x's) overlayed with (j)n{r) 
obtained by shooting equations 20 and 21 (solid curves). The Fourier-decomposed PDE data overlays the shooting solution 
everywhere. 
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FIG. 10: Power spectra of the core amplitude, (j){t, 0), for the oscillons barely above and below the ro ~ 2.335 resonance. The 
power measured in each frequency regime slowly diminishes as the oscillon radiates away much of its energy until approximately 
t = 1100 where the oscillon enters a non-radiative state and all the components of the power spectrum become constant. 
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FIG. 11: Plot of 0) for ro = 7.25 displaying extremely nonlinear and unpredictable behavior during the "bouncing" 
phase (for t < 60), after which the field settles into a typical oscillon evolution. Once in the oscillon regime, the period is 
approximately T ~ 4.6, and the first two modulations of the field can be seen (envelope maxima at f « 105 and t ~ 200). 
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FIG. 12: Plot of oscillon lifetime versus initial radius of bubble for 4.22 < ro < 9. Although there seem to be no oscillons in 
the range 4.6 ^ ro ^ 6, it is clear that oscillons and resonances do exist for higher initial bubble radii, ro ^ 6.5. 
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FIG. 13: Fundamental field (/!>(r) in the freeze-out region at i = 0, 57, 75, and 270. The characteristic speeds of the radiation 
as r —> r-c (here, Tc = 56) and the wavelength of the radiation is blue-shifted to the lattice Nyquist limit, 2Ar. The 
Kreiss-Oliger dissipation explicitly added to the finite difference equations subsequently "quenches" the field. 
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FIG. 14: Convergence factor, C/ = ||04h — 4>2h\\2/\\4>2h — 4>h\\2, for the field 4> composed from the solutfon at three different 
discretizations (value of 4 indicates 2nd order convergence). The I2 norm || ■ ■ • II2 is defined by \\v\\2 ~ i^N~^ "Yln-i '^O 
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FIG. 15: Plot of energy contained in oscillon (dashed line), energy radiated (dotted line), and total energy (solid line). The 
total energy of the system is a constant of motion and is numerically conserved to within a few tenths of a percent. The 
energy contained within the bubble drops rapidly during the initial radiative phase and plateaus around E ~ 43m/ A during 
the quasi-stable "oscillon" phase. 
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FIG. 16: Plot comparing the OBC (solid line) and MIB (dashed line) solutions to an "ideal" solution. The OBC solution is 
obtained using a massless outgoing boundary condition, the MIB solution is obtained by solving the system in spherical MIB 
coordinates, and the ideal solution is obtained by evolving the solution in standard {r,t) coordinates on a grid large enough 
to ensure no reflection off the outer boundary. The error estimates are obtained from the ^2-norm of the difference between 
the trial solutions (OBC or MIB) and the ideal solution, 110 — ^idoai 1 1 2 ■ Contamination of the OBC solution is observed at two 
crossing times, t « 120, where the error estimate increases by over three orders of magnitude. The MIB solution error grows 
slowly and steadily as expected when solving a continuum equation with two different finite difference techniques. 
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FIG. 17: 0) versus time for the large-grid solution (triangles), MIB solution (solid curves), and OBC solution (dashed 

curves). The solutions all agree before 2tcrossing, but the OBC solution begins to drift away from the ideal solution after 
2tcrossing. The crror in the OBC solution is due to radiation that is reflected off of the outer boundary (hence needing two 
crossing times to return to r = to contaminate the oscillon). All pictures span the same area, A<j!> — 0.075 by At = 0.5. 



